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Abstract 

A computer adapted fluctuation formula for the calculation of the wavevec- 
tor- and frequency-dependent dielectric permittivity for interaction site mod- 
els of polar fluids within the Ewald summation technique is proposed and 
applied to molecular dynamics simulations of the TIP4P water. The formula 
is analyzed and optimal parameters of the Ewald method are identified. A 
comparison of the obtained results with those evaluated within the reaction 
field approach is made. 
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1 Motivation 



In order to achieve a macroscopic behaviour for investigated quantities in com- 
puter experiment based on the observation of finite systems, it is necessary to reduce 
the influence of surface effects to a minimum. This is especially important for polar 
systems with the long-range nature of interactions. Excluding the surface effects in 
simulations can be performed within either the reaction field (RF) [1-5] or Ewald 
summation [6-10] techniques. Now an equivalence between these techniques has 
been established for models of point dipoles and proper calculations can be made 
within either method [11, 12]. The explicit consideration of a finite-size medium 
lead to computer adapted fluctuation formulas [11-17] which allow one to calculate 
boundary free values for the dielectric constant on the basis of dipole moment fluc- 
tuations obtained in simulations. These formulas differ considerably with respect to 
those known from the theory of macroscopic systems even if the Ewald method is 
used [11]. Details of the summation must be taken into account explicitly in order 
to obtain correct values for the bulk dielectric constant. 

Previously [18-20], the standard RF of point dipoles (PDRF) [3] has been applied 
to investigate more realistic, interaction site models (ISMs) [21] of polar fluids. 
The PDRF, however, being exact for point dipole models, may not be necessarily 
applicably to interpret simulation results for arbitrary systems [5]. Recently, it has 
been shown by actual calculations for a MCY water model that uncertainties for the 
dielectric quantities are significant if the PDRF is used in computer simulations of 
ISMs and an alternative scheme, the interaction site reaction field (ISRF) geometry, 
has been proposed [22]. At the same time, there is not such an approach concerning 
the entire wavevector and frequency dependence for the dielectric permittivity of 
ISMs within the Ewald geometry. The main attention of previous simulations [23- 
31] was directed to study the dielectric properties in the static limit or at zero 
and small wavevector values. Moreover, the macroscopic fluctuation formulas have 
been used in the simulation results without taking into account details of the Ewald 
summation. 

In the present paper we apply the Ewald technique for treating Coulomb inter- 
actions in ISMs. The paper is organized as follows. A fluctuation formula suitable 
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for the calculation of the wavevector- and frequency-dependent dielectric constant 
is derived in Sec. 2 and optimal values of the Ewald parameters are determined 
there. The results of molecular dynamics simulations of the TIP4P water for time 
correlation functions related to dielectric polarization are presented in Sec. 3. These 
results are compared with those computed within the ISRF geometry. Concluding 
remarks are given in Sec. 4. 



2 Ewald summation for ISMs 

Consider a polar fluid with N molecules composed of M interaction sites which 
are confined in a volume V. The microscopic electric field created by the molecules 
at point r and time t can be presented as E(r,t) = / D(r — r')Q(r', t)dr', where 
Q(r,t) = Y J f=iY J ^=iQa^{ r — r i(t)) is the microscopic charge density, r?(t) and q a 
denote the position and charge, respectively, of site a within the molecule % and 
D{p) = —V 1/p is the operator of the Coulomb interactions. 

Obviously, the field E(r,t) for infinite systems (N,V — > oo) can not be repro- 
duced exactly in computer experiment which deals with a finite, as a rule, cubic 
volume V = L 3 , where L is the length of the simulation box edge. However, using 
the lattice summation, a macroscopic behaviour can be achieved considering the 
interactions between sites within the basic cell as well as an infinite lattice of its 
periodic images (the periodic boundary convention). This can be interpreted as 
an effective interaction which involves only the sites in the basic cell and charac- 
terized by a modified operator T>(p) = J2 n D(p + nL), where the summation is 
extended over all vectors n with integer components. It is more convenient to rep- 
resent the lattice sum in a form, proposed by Ewald and Kornfeld (EK) [6], namely, 
V(p) = Di(p) + D 2 (p), where 

2r] 

u{p-\-nij)<.eric{ri\p-\-n±j\ + 

0<\n\<Af 

is a sum in real coordinate space, while 



Di{p)= D(p + nL)[ei{c(r]\p^nL\ + ^\p + nL\e^(-rj 2 \p + nL\ 2 )] 

0<\n\<M V 71 " 



D2(P) = ^ E D(k)eM-k 2 /±V 2 + ik.p) (2) 

0<|fc|<fc max 



corresponds to summation over wavevectors k = 2nn/L of the reciprocal lattice 
space and D{k) = j dre~ lk ' p D(p) = —Anik/k 2 is the spatial Fourier transform 
of D(p)- For the idealized summations (A/" — > oo, k ma , x — > oo), the total sum of 
(1) and (2) is independent on the parameter 77. The main advantage of the EK 
representation lies in the fact that values for 77 can be found in such a way that the 
both sums, Di and D 2 , converge very quickly and may be truncated after a finite 
number of terms. If the parameter 77 is chosen sufficiently large, we can restrict 
ourselves to a single term (A/" = 0) in the real space sum, corresponding to the 
basic cell to which toroidal boundary conditions are applied, and, additionally, to 
the spherical truncation \p\ < R, where R < L/2. 

In such a case, taking the Fourier transforms of (1) and (2), after some algebra 
one obtains -Di(fc) = —A7riD 1 (k)k/k 2 and D 2 {k) = —A-K\D 2 {k)k/k 2 , where 

£>!(*) = j k 3l (kp)(eM V p) + ^lpexp(-7 ? V))dp , (3) 

D 2 (k) = exp(—k 2 /4r] 2 ) if < k < k max and D 2 {k) = otherwise and j^z) = 
sin(,2)/,2 2 — cos(z)/z denotes the spherical Bessel function of first order. Then the 
Fourier transform of the electric field is 

E{k, t) = (Di(fe) + D 2 (kj)Q{k, t) = -AnP L (k, t)D(k) , (4) 

where Q(k,t) = Ej M ^e" ifc ' r?W , P L (JM) = %Q(k,t) is the longitudinal compo- 
nent of the microscopic operator P of polarization density (V-P(r,t) = — Q(r,t)) 
and D(k) = D^ + D^k). 

Let us apply an external electric field E Q {k,uj) to the system under consider- 
ation. The longitudinal, wavevector- and frequency- dependent dielectric constant 
is defined via the material relation 4irPi,(k,uj) = (e L (k,u) — ljE L (k,uj), where 
P L (fc,c<j) = ^P L (fc,w)^ and Ei,(k,u) = (Jek'E Q (k,u) + E(k,uj)^ are macroscopic 
values for longitudinal components of the polarization and total field, ( ) denotes 
statistical averaging at the presence of the external field and the time Fourier trans- 

/oo 
dte~ iu}t J?(k,t) has been used for the functions P L (k,t), 

E(k,t) and k = k/k. Perturbation theory of the first order with respect to E Q 

4 



yields P L (k,u) = - v ^JJ > dte-^(P L (k,OyP L (-k,t)) Q kk'E (k,u), where 
( ) denotes equilibrium averaging at the absence of the external field, and k-Q and 
T are the Boltzmann's constant and the temperature of the system, respectively. 
Then, eliminating E Q (k,u) from the presented above expressions, we obtain the 
desired fluctuation formula 



s&g - 1 = ^(-Wfh = 9y ^ ( _ {K t)) . (5) 



Here G^k.t) = (^Pi J (k,0)-P^(—k,t)^^Nfi 2 is the longitudinal component of the 
wavevector-dependent dynamical Kirkwood factor for the finite system, jj, = I/Ltj 

denotes the permanent magnitude of the molecule's dipole moment /it, = 1a r ii 

i f°° 
y = AnNfi 2 /9Vk B T and 3? iui (...) = J ... e~ lw 'dt is the Laplace transform. The 

right-hand side of Eq. (5) corresponds to the well-known fluctuation formula for 

infinite systems, where g h (k,t) = liniTv^oo G^(k, t) is the infinite-system Kirkwood 

factor. 

The computer adapted formula (5) reduces to the formula for infinite systems 
if the function D{k) = 1. It can be shown easily that for nonzero wavevectors the 
function D{k) — > 1 if /c max — > oo, additionally provided r\ — > oo at M = 0. For 
k = the pattern is different because of finiteness of L and D(0) — as in the ISRF 
geometry [22]. However, in the case of an actual summation, when k max takes finite 
values, the factor D(k) can noticeably differ from unity. Therefore, the finite sample 
behaves like a macroscopic system if the function D{k) is very close to unity and 
this condition can be verified now quantitatively. Moreover, this explicit result may 
serve as an initial point for a more fruitful discussion about the Ewald method itself. 
Let A = maxfc|l — D(k)\ be a maximal deviation of D(k) from unity in the whole 
interval of acceptable nonzero wavevector values for a chosen pair of parameters r\ 
and /c max - Then an optimal value for r\ can be determined as that providing a global 
minimum for the function A(i], fc max ) at a given k max . 

According to Eq. (5), the obtained in simulations Kirkwood factor Gx differs 
from its genuine value g L with the relative precision of x — 9yA. The function 
x(Vi n max) is shown in Fig. 1 as depending on r\ at fixed values of n max = k mauX L/2n. 
It has been calculated for the case of R = L/2 and y = 5.47 that corresponds to 



the thermodynamics point p = mN/V = 1 g/cm 3 , T = 293 K of the TIP4P model, 
where m is the mass of water molecule. As we can see from the figure, the precision 
of calculations of dielectric quantities in computer experiment depends on Ewald 
parameters in a characteristic way. We indicate the existence of the sharp minimum 
of x at an arbitrary value of n max . The curves of Fig. 1 can be useful to estimate 
the possibility of a given simulation result to reproduce directly the macroscopic 
dielectric behaviour of an IS system in an arbitrary thermodynamics state, because 
then the function x' — XU'/u is simply rescaled, using the actual value of y'. From 
the last equality it follows that the precision of calculations is better for systems 
with lower particle densities N/V, molecular polarities \i and higher temperatures 
T. It is obvious also that minimums of the functions A and x with respect to Ewald 
parameters coincide between themselves. 

The optimal pairs of values for rj and n max at R = L/2 as well as the cor- 
responding values of the functions A and x are selected in Table 1. Choosing a 
criterion x < 1%, we may ask that the formula for infinite systems might be applied 
(at k 7^ 0) and the influence of summation details can be neglected in this case for 
which Gi,(k, t) and g L (k, t) are indistinguishable. It can be seen easily from the table 
that values of n max > 4 satisfy this criterion if the parameter 77 is chosen optimally. 
The parameters n max = 5, r\L = 5.76 and R = L/2 are usually exploited in simula- 
tions [10]. For these values the relative precision is x — 0.22%. However, choosing 
the optimal value r\L = 5.929 at n max = 5 instead of r\L = 5.76, we can reduce the 
uncertainty up to x — 0.13%. 

In the presented above consideration, the cut-off radius R has been putted to 
be half the basic cell length. Nevertheless, increasing n max , the same precision of 
summation can be achieved also at smaller values of R. Let i] and n max correspond 
to the optimal parameters at R = L/2. And now we choose a smaller value of the 
cut-off radius in the form R' = R/l, where I > 1. Taking into account the fact that 
maximum deviations of D(k) from unity are always observed at k — 27r(n max + l)/L, 
it is easy to show that the same value of A(?7, n max ) can be obtained also at rj' = Irj 
and n^ ax = Z(n max + 1) — 1. For example, putting n max = 5 and r/L = 5.929 at 
R = L/2, we then obtain for R = L/A{1 = 2) the following results: i]'L = 11.858 and 
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n max — 11- Choosing smaller values of R can be more convenient if the summation 
in fc-space takes less computation time in an actual programme than the summation 
in real space. Indeed, let t\ and t 2 are the computation times in real and fc-space 
(ti > t 2 ), respectively, at given values of R and n max . It is obvious that t 1 ~ R 3 and 
h ~ (rw + l) 3 - Then using new values n' max = Z(n max + 1) — 1 and R' = R/l and 
minimizing the total computation time t' = t[ + 1' 2 with respect to /, one obtains 
I = \fhjt~2- Therefore, in such a way we can provide even a time optimization of 
the programme without any loss of the precision. 

3 Numerical results. Comparing the Ewald and 
reaction field methods 

The study of dielectric properties by computer experiment is still a major chal- 
lenge, given that the calculations are very sensitive to long-range interactions and 
because the polarization of polar fluids is a collective effect, so that long trajectories 
are required in order to obtain adequate statistical accuracy. For this reason, un- 
til now, the dynamical polarization of ISMs has been investigated at zero or small 
wavevector values only [18, 19, 23, 24, 27, 29, 31]. As far as we know, there are 
no computer experiment data on the entire wavevector dependence of dynamical 
dielectric quantities for such systems. 

Our molecular dynamics simulations were carried out for the TIP4P model [32] 
in the microcanonical ensemble at a density of p = 1 g/cm 3 and at a temperature 
of T = 293 K. We have performed two runs corresponding the Ewald and ISRF [22] 
geometries, respectively. In the both runs N = 256 molecules were considered in the 
cubic sample V = L 3 to which toroidal boundary conditions were applied (A/" = 0) 
and the interaction cut-off radius was half the basic cell length, R = L/2 = 9.856A. 
The simulations were started from a well equilibrated configuration for positions 
of sites, obtained by Monte Carlo simulations. Initial velocities of molecules were 
generated at random with the Maxwell distribution. The equations of motion were 
integrated with a time step of At = 2 fs on the basis of a matrix method [33] using 
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the Verlet algorithm in velocity form. The system was allowed to achieve equilibrium 
for 50 000 time steps. The equilibrium state was observed during 500 000 At = 1 ns 
and each 10th time step was chosen to compute equilibrium averages. Translational 
and angular velocities of molecules were slightly rescaled after every 500 time steps 
in order to conserve the total energy of the system, so that the relative total energy 
fluctuations did not exceed 0.01% over the whole runs. 

The dynamical Kirkwood factor was evaluated in the time interval of lOOOAt = 2 
ps and in a very large wavenumber region, namely, at k — [0, 1, ... , 300] A; min , where 
^min = 2n/L = 0.319A \ Considering the system during such a rather long period 
of time allows us to achieve statistical accuracy for the investigated quantities of 
order 1%. The optimal parameters i]L = 5.929 and n max = 5 have been used 
in the Ewald summation of Coulomb forces. The computational times on IBM PC 
AT486DX4 100 MHz to evaluate dynamics of the system in our Fortran programmes 
were 2.2 s and 1.2 s per step in the cases of Ewald and ISRF geometries, respectively. 

Within the Ewald geometry the Coulomb part q a Q b /\Pij\ of the intersite potential 
is replaced by 

k = *4 ( r - w) E^«»(*-p?)}=^(W)+^(pff)- 

' \Pij\ V \k\>0 fc J 

(6) 



Here, p ab , = r" — r b designates the distance between sites belonging the basic cell 
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(rf,r b G V), J}"!- = r" — r b , where, according to the toroidal boundary conditions, 
r b = r b +pL (p = (p x ,Py,Pz)', Px,P y ,Pz = 0, ±1) is the position of the nearest image 
of r b with respect to rf , and 9 denotes the Heviside function, i.e., 6(p) = 1 if p > 
and 9(p) = otherwise. The function 6 indicates about the spherical site-site 
truncation in the real coordinate space. The force acting on the a-th charged site 
of molecule i due to the interaction with the 6-th charge of molecule j is F® b = 
—d(p ab / dr°" or in a more explicit form 



erfc^l^D + ^lp^le- 2 ^! 2 



4-7T ^ ma ' x \c ^ 

+ T7 E V ^ Sin(fc.pg') EE q a q b (Dm-\) + D 2 {p<* 



(7) 



—ab 

We note that the 5-like part q a q b ! ^-S(R— \p^\)eric(rfR) of the force is not included 
in (7) for the reason that the complementary error function vanishes at \pfj\ = R 
for sufficiently large values of r\. In particular, in our case erfc^i?) = erfc(?iL/2) = 
0.0000276 <C 1. Then the potential (6) can be considered as a continuous and 
continuously differentiable one and the drift of the total energy of the system, as- 
sociated with the passage of sites through the surface of the truncation sphere, can 
be neglected. 

In an actual molecular dynamics programme the current potential energy of 
the system, U = \ J2i^j Z)£f& an d the total force acting on the a-th site of 
molecule % due to interactions with all the rest of sites belonging other molecules, 
F1 = Ef=i ££Li F$, can be calculated as follows 



N M 1 n max p-^TZT N M 

EE^(Ip?I) + ^ E -^r^(n)f(-n)-EE^(p") 

i^j a,b \n\>0 i=l a,b 



, (8) 



N M M 



9(7 ™ max Tl TT 2 n 2 

^=^E E^l(l7^|) + ^E 7^ 

3=1 6=1 ^ |ri|>0 71 6=1 

(9) 



where the self electrostatic energy, u — | Yh=i YJa^b ^ >< f i -> an d the self forces, / ® = 
J2bLiFiii have been excluded from (8) and (9), respectively, because the intramolec- 
ular forces do not contribute into molecular translational accelerations and torques. 
The auxiliary function 

N,M 

<3f{n) = E Vae~ 2nin< /L = Re^(n) + i Im^(n) (10) 



is introduced in order to reduce the total number of numerical operations in k- 
space from of order (A r M) 2 n max to iVMn max that is very important for simulat- 
ing large systems. This number can be reduced approximately twice yet, using 
invariance of the subsume expressions with respect to the inverse transformation 
n — > —n. Finally, taking into account that the real part of W (n) is an even 
function of n and the imaginary part is an odd one, we obtain that only the real 
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part Reie- 2 ™- r "/ L ^(-n) = sm(2irn-r° / L)Re& (n) + cos(27rn-r?/L)Iin0 r (n) give 
nonzero contributes into (9) and (n)<3f (-n) = (Re^(n)) 2 + (Im^(n)) 2 . 

In the RF geometry, real particles of the infinite system, which are located outside 
the sphere of finite radius R around a reference particle belonging the basic cell, are 
replaced by an infinite, as a rule, conducting continuum. There are two versions of 
the RF geometry. In the PDRF approach, molecules are considered as point dipole 
particles and the intermolecular potential is of the form [18, 19]: 



where Vi is the centre of mass of the i-th molecule and the molecular cut-off is 
performed. In the exact ISRF method [22] the spatial distribution of charges within 
the molecule is taken into account explicitly at constructing the reaction field. As a 
result, the potential (11) transforms into 



where the first term in the right-hand side of (12) describes the usual Coulomb field, 
whereas the rest of terms corresponds to the reaction field in the IS description. 

It can be shown easily that the potential (12) is reduced to (11) in one case only, 
namely, when d/R — > 0, where d = 2max a |r" — Ti\ denotes the diameter of the 
molecule. In this case, the positions for sites and centres of mass are undistinguished 
within the same molecule. For finite samples of IS molecules we have d/R ^ and, 
therefore, the PDRF potential (11) may affect on a true macroscopic behaviour of 
the system considerably. Moreover, the ISRF method has yet a minor advantage 
over the PDRF scheme that the potential of interaction (12) is continuous and 
continuously differentiable. It is worth to mention also that in the RF geometry the 
dielectric permittivity is computed using the fluctuation formula (5) with the formal 
substitution D(k) -> D RF {k) = 1 - Sj^kR) / \kR) [22]. 

The wavevect or- dependent static Kirkwood factor, G^(k) = Gl(A;,0), and sam- 
ples of the normalized dynamical Kirkwood factor, <3> L (/c,t) = G L (k,t)/G L (k), cal- 
culated in the simulations within the Ewald and ISRF geometries, are shown in 




(11) 




(12) 
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Figs. 2, 3 by the circles and dashed curve, respectively. Since in the ISRF geometry 
the function -Drf(&) differs from unity considerably, to evaluate the infinite system 
Kirkwood factor g L (k,t) the performance of the self-consistent transformation (5) 
is necessary. This result is plotted by the solid curve. At the same time, within 
the Ewald geometry the function D{k) is very close to unity at the given optimal 
parameters of summation (see Table 1), so that the infinite system Kirkwood factor 
is equivalent to that, obtained directly in the simulations, i.e, g h (k, t) = G L (k, t) (ex- 
cepting the case k — 0). As we can see from the figures, the agreement between the 
two sets of data for the infinite-system functions, corresponding to the Ewald and 
ISRF geometries, is quite good. The slight difference (within a few per cent) at large 
times can be explained by an approximate character of the integration appearing 
for the ISRF geometry at performing the inverse Laplace transform of (5). 

For the purpose of comparison, the infinite-system Kirkwood factor g L (k) corre- 
sponding to the PDRF geometry is also included in the Fig. 2 (the dotted curve). 
Deviations of values for g L (k) obtained using the PDRF potential from those eval- 
uated in the Ewald and ISRF geometries are of order 20%. They are well exhibited 
at intermediate values of wavevectors. Such a situation can be explained by the 
fact that the PDRF geometry does not take into account the spatial distribution of 
charges within the molecule and, thus, the precision of calculations for wavevector- 
dependent dielectric quantities at k ~ 2rr/d ~ 3.4A 1 can not exceed d/R ~ 20%, 
where d = 1.837A for the TIP4P water molecule. And only for great wavevec- 
tor values (k > 6A 1 ), where the influence of boundary conditions is negligible 
(.Drf(^) — ► 1), all the three geometries become completely equivalent. 

4 Conclusion 

Explicitly considering details of the Ewald summation to treat Coulomb interac- 
tions, the fluctuation formula for the computation of the dielectric permittivity in IS 
models of polar fluids has been rigorously derived. Using this formula, it has been 
corroborated by actual molecular dynamics calculations that the Ewald and ISRF 
methods can be applied with equal successes to investigate the dielectric constant 
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of ISMs in computer experiment. The Ewald geometry, however, at a specific choice 
for parameters of the summation, may reproduce the macroscopic behaviour for 
dielectric quantities directly in simulations without any additional transformations. 

Since the calculation of the wavevector- and frequency- dependent dielectric per- 
mittivity in simulations for ISMs is practical now in principle, we believe that this 
fact will stimulate further research of such systems in theory, computer and pure 
experiment. 
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Table 1. Optimal parameters of the Ewald summation for ISMs at R — L/2 



riL 




A(%) 


x(%) 


3.301 


1 


1.632E+00 


80.29 


3.874 


2 


5.523E-01 


27.17 


4.791 


3 


6.684E-02 


3.288 


5.209 


4 


2.212E-02 


1.088 


5.929 


5 


2.690E-03 


0.1324 


6.276 


6 


8.978E-04 


0.0442 


6.887 


7 


1.094E-04 


0.0054 


7.251 


8 


3.602E-05 


0.0018 



Figure captions 

Fig. 1. The precision of reproducing bulk dielectric quantities in computer ex- 
periment for ISMs as depending on parameters of the Ewald geometry. 

Fig. 2. The static wave vector- dependent Kirkwood factor of the TIP4P water. 
The obtained result in the Ewald geometry is presented by the full circles. The 
dashed and solid curves correspond to the finite and infinite systems in the ISRF 
geometry. The PDRF infinite-system Kirkwood factor is shown as the dotted curve. 

Fig. 3. The normalized dynamical Kirkwood factor of the TIP4P water. Nota- 
tions as for fig. 2. 
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